{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Getting the necessary data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You just need to do this only once"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "!rm -f genotypes.vcf.gz 2>/dev/null\n",
    "!tabix -fh ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/release/20130502/supporting/vcf_with_sample_level_annotation/ALL.chr22.phase3_shapeit2_mvncall_integrated_v5_extra_anno.20130502.genotypes.vcf.gz 22:1-17000000|bgzip -c > genotypes.vcf.gz\n",
    "!tabix -p vcf genotypes.vcf.gz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/tiago/anaconda3/lib/python3.5/importlib/_bootstrap.py:222: RuntimeWarning: numpy.dtype size changed, may indicate binary incompatibility. Expected 96, got 88\n",
      "  return f(*args, **kwds)\n",
      "/home/tiago/anaconda3/lib/python3.5/importlib/_bootstrap.py:222: RuntimeWarning: numpy.dtype size changed, may indicate binary incompatibility. Expected 96, got 88\n",
      "  return f(*args, **kwds)\n"
     ]
    }
   ],
   "source": [
    "from collections import defaultdict\n",
    "\n",
    "%matplotlib inline\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import vcf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Variant Level information\n",
      "CIEND\n",
      "CIPOS\n",
      "CS\n",
      "END\n",
      "IMPRECISE\n",
      "MC\n",
      "MEINFO\n",
      "MEND\n",
      "MLEN\n",
      "MSTART\n",
      "SVLEN\n",
      "SVTYPE\n",
      "TSD\n",
      "AC\n",
      "AF\n",
      "NS\n",
      "AN\n",
      "ASN_AF\n",
      "EUR_AF\n",
      "AFR_AF\n",
      "AMR_AF\n",
      "SAN_AF\n",
      "DP\n",
      "Sample Level information\n",
      "GT\n",
      "DP\n"
     ]
    }
   ],
   "source": [
    "v = vcf.Reader(filename='genotypes.vcf.gz')\n",
    "\n",
    "print('Variant Level information')\n",
    "infos = v.infos\n",
    "for info in infos:\n",
    "    print(info)\n",
    "\n",
    "print('Sample Level information')\n",
    "fmts = v.formats\n",
    "for fmt in fmts:\n",
    "    print(fmt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "22 16050075 None A [G] 100 []\n",
      "{'NS': 2504, 'SAN_AF': [0.0], 'AN': 5008, 'AF': [0.000199681], 'AMR_AF': [0.0], 'EAS_AF': [''], 'DP': [8012], 'EUR_AF': [0.0], 'ASN_AF': [0.0], 'SAS_AF': ['0.0010'], 'AC': [1], 'AFR_AF': [0.0]}\n",
      "GT:DP\n",
      "2504\n",
      "True ['0', '0'] False False True\n",
      "1\n"
     ]
    }
   ],
   "source": [
    "v = vcf.Reader(filename='genotypes.vcf.gz')\n",
    "rec = next(v)\n",
    "print(rec.CHROM, rec.POS, rec.ID, rec.REF, rec.ALT, rec.QUAL, rec.FILTER)\n",
    "print(rec.INFO)\n",
    "print(rec.FORMAT)\n",
    "samples = rec.samples\n",
    "print(len(samples))\n",
    "sample = samples[0]\n",
    "print(sample.called, sample.gt_alleles, sample.is_het, sample.is_variant, sample.phased)\n",
    "print(int(sample['DP']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "defaultdict(<class 'int'>, {('sv', 'DEL'): 6, ('indel', 'del'): 273, ('snp', 'ts'): 10054, ('indel', 'ins'): 127, ('snp', 'tv'): 5917, ('sv', 'CNV'): 2, ('snp', 'unknown'): 79, ('sv', 'SVA'): 1, ('indel', 'unknown'): 13})\n",
      "defaultdict(<class 'int'>, {1: 15971, 2: 79})\n"
     ]
    }
   ],
   "source": [
    "f = vcf.Reader(filename='genotypes.vcf.gz')\n",
    "\n",
    "my_type = defaultdict(int)\n",
    "num_alts = defaultdict(int)\n",
    "\n",
    "for rec in f:\n",
    "    my_type[rec.var_type, rec.var_subtype] += 1\n",
    "    if rec.is_snp:\n",
    "        num_alts[len(rec.ALT)] += 1\n",
    "print(my_type)\n",
    "print(num_alts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = vcf.Reader(filename='genotypes.vcf.gz')\n",
    "\n",
    "sample_dp = defaultdict(int)\n",
    "for rec in f:\n",
    "    if not rec.is_snp or len(rec.ALT) != 1:\n",
    "        continue\n",
    "    for sample in rec.samples:\n",
    "        dp = sample['DP']\n",
    "        if dp is None:\n",
    "            dp = 0\n",
    "        dp = int(dp)\n",
    "        sample_dp[dp] += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0x7fa4cfc124a8>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA8EAAAIMCAYAAADVSa7eAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3Xmc1lXB///XYRNcQIRBEVBccCdRySUtt0TKCiMiTQcUDMslSyvRFvt1W2mLlpXeqLmngma5Y2bAqCGKzjjgjjuug+AGynp+f3w+fO9RB2ZhZs61vJ6Pxzyu4Vyf6/q8se770dtzPueEGCOSJEmSJJWDDqkDSJIkSZLUXizBkiRJkqSyYQmWJEmSJJUNS7AkSZIkqWxYgiVJkiRJZcMSLEmSJEkqG5ZgSZIkSVLZsARLkiRJksqGJViSJEmSVDYswZIkSZKkstEpdYD20rt37zhw4MDUMVQEnqtbDMDWFRskTiJJkiSpqR5++OEFMcaKxq4rmxI8cOBAZs+enTqGisA3Js0EYPLx+yROIkmSJKmpQggvNuU6l0NLkiRJksqGJViSJEmSVDYswZIkSZKksmEJliRJkiSVDUuwJEmSJKlsWIIlSZIkSWXDEixJkiRJKhuWYEmSJElS2bAES5IkSZLKhiVYkiRJklQ2LMGSJEmSpLJhCZYkSZIklY1GS3AIoWsI4cEQwqMhhMdCCP9fPn5FCOH5EEJN/jMkHw8hhAtCCPNCCLUhhN3rfdfYEMIz+c/YeuN7hBDm5J+5IIQQ8vFNQgh359ffHULo2dg9JEmSJElak6bMBC8FDoox7goMAYaHEPbO3/thjHFI/lOTj30BGJT/TAAugqzQAmcBewF7AmetLrX5NRPqfW54Pj4RuCfGOAi4J//zGu8hSZIkSdLaNFqCY+b9/I+d85+4lo+MAK7KP/cAsHEIoS9wKHB3jHFhjHERcDdZoe4LdI8xzowxRuAq4PB633Vl/vuVHxtv6B6SJEmSJK1Rk54JDiF0DCHUAG+SFdlZ+Vu/zJcjnx9CWC8f6we8XO/j8/OxtY3Pb2AcYNMY42sA+WufRu7x8dwTQgizQwiz6+rqmvJXlSRJkiSVsCaV4BjjyhjjEKA/sGcIYRfgDGAH4NPAJsDp+eWhoa9owfjaNOkzMcaLY4xDY4xDKyoqGvlKSZIkSVKpa9bu0DHGt4HpwPAY42v5cuSlwOVkz/lCNis7oN7H+gOvNjLev4FxgDdWL3POX99s5B6SJEmSJK1RU3aHrgghbJz/3g34PPBkvXIayJ7VnZt/5BZgTL6D897AO/lS5ruAYSGEnvmGWMOAu/L33gsh7J1/1xjg5nrftXoX6bEfG2/oHpIkSZIkrVGnJlzTF7gyhNCRrDRPiTHeFkL4Twihgmxpcg3w7fz6O4AvAvOAJcCxADHGhSGE/wEeyq/7RYxxYf77d4ArgG7AnfkPwDnAlBDCeOAl4Otru4e0Tl55BW68EZ7oDDHCW9tBr16pU0mSJElqRSHbkLn0DR06NM6ePTt1DBWa11/Piu/kyXDffQB847gL4IMlTJ73D/j3v2HDDROHlCRJktSYEMLDMcahjV3XrGeCpZLw5pvwv/8LBx4Im28OJ58Mb78Nv/gFPPkkDB0KO+4EDz0Eo0bBsmWpE0uSJElqJZZglYe33oJLLoFDDoG+feE738lmgX/2M5g7F+bMgZ/+FLbfPru+d+/s+rvugmOOgVWrksaXJEmS1Dqa8kywVJwWLYJ//jNb6vzvf8PKlbDttnDGGfCNb8Auu0Bo6LSt3LhxUFcHEydmpfiPf1z79ZIkSZIKniVYpefGG+GKK+Bf/4Lly2GrreAHP8iK75AhzSuyP/pRtnz6vPOgTx/4yU/aLLYkSZKktmcJVmm56y74+tdhiy3glFNg9OjsGd+WzuCGAL/9LSxYkC2X7t0bvv3txj8nSZIkqSBZglVazjkH+vWDZ56BLl1a5zs7dIBLL82eKz7hhKwIjxrVOt8tSZIkqV25MZZKx6xZMH06nHpq6xXg1Tp3hilT4DOfgaOOgnvuad3vlyRJktQuLMEqHeeeCz17wre+1Tbfv/76cOut2Q7Shx8OnjstSZIkFR1LsErDE0/AP/4BJ50EG23Udvfp2ROmTs2WRH/hC/D00213L0mSJEmtzhKs0vDb30K3bnDyyW1/r803z3aeDgGGDYNXXmn7e0qSJElqFZZgFb/58+Gaa2D8eKioaJ97DhqUzQgvXAiHHpq9SpIkSSp4lmAVv/PPh1Wr4LTT2ve+u+8ON9+c7UT9pS/BkiXte39JkiRJzWYJVnFbuBAmTYIjj4SBA9v//gceCNddl+1MPWoULF/e/hkkSZIkNZklWMXtwgth8WL40Y/SZRg5Ev73f+HOO2HcuGxWWpIkSVJB6pQ6gNRiS5bAH/8Ihx0GgwenzfKtb0FdHfz4x9nO0eedl22cJUmSJKmgWIJVvC67DBYsgIkTUyfJnHEGvPkm/OEP0KdP9mdJkiRJBcUSrOK0fDn87nfwmc/AfvulTpMJIZsBXrAAzjwz20F61KjUqSRJkiTV4zPBKk5TpsCLLxbOLPBqHTrA5ZfDrrvCT34CK1emTiRJkiSpHkuwik+McO65sPPO2fPAhaZz5+zZ4Keegn/8I3UaSZIkSfVYglV87rwT5syB00/PZl4L0ciRsN128KtfZaVdkiRJUkEo0AYhrcU558AWW8ARR6ROsmYdO2ZLtaurYerU1GkkSZIk5SzBKi733w/33gunnZYtOy5kRx+dlfVf/Sp1EkmSJEk5S7CKy7nnQq9eMH586iSN69wZfvhDuO8+qKpKnUaSJEkSlmAVk8ceg1tvhZNPhg02SJ2macaPz84MdjZYkiRJKgiWYBWP3/wG1l8fTjopdZKm69YNTj0V7roLHn44dRpJkiSp7FmCVRxeegmuvRYmTMiWQxeT73wHevRwNliSJEkqAJZgFYfzzsteTz01bY6W6N49W8J9003w+OOp00iSJEllzRKswrdgAVxyCRx1FAwYkDpNy5xySraU+5xzUieRJEmSypolWIXvz3+GJUvgRz9KnaTleveG44/PlnQ//3zqNJIkSVLZsgSrsC1eDH/6E4wYATvtlDrNujntNOjYEX7729RJJEmSpLJlCVZhu/RSWLgQTj89dZJ1168fHHMMXHYZvPZa6jSSJElSWbIEq3AtXw6//z187nOwzz6p07SOH/0o+3ut3uhLkiRJUruyBKtwXXcdvPwyTJyYOknr2WYbOPJIuOiibIZbkiRJUruyBKswrVoF554Ln/oUDB+eOk3rmjgxe9b5ggtSJ5EkSZLKjiVYhem227IzdU8/HUJInaZ17bJLttHXBRfAe++lTiNJkiSVFUuwCk+M2Xm6AwfC6NGp07SNM8+ERYtg0qTUSSRJkqSyYglW4bnvPpg5E374Q+jUKXWatrHnnvD5z2cbf334Yeo0kiRJUtmwBKvwnHMOVFTAscemTtK2zjwTXn8dLr88dRJJkiSpbFiCVVhqa+GOO+CUU6Bbt9Rp2tYBB2RHP517bnZskiRJkqQ2ZwlWYTn/fNhwQzjhhNRJ2l4I2Wzwiy9mx0FJkiRJanOWYBWOGLNZ4BEjoGfP1Gnax2GHZcdA/frX2bFQkiRJktqUJViF44kn4M034cADUydpP6tng598Ev7xj9RpJEmSpJJnCVbhmD49ez3ggJQp2t+oUbDttvCrX2Wz4ZIkSZLajCVYhWP6dBgwALbeOnWS9tWxI0ycCI88Av/6V+o0kiRJUkmzBKswxJiV4AMOyJYIl5vKSujfH375y9RJJEmSpJJmCVZhePxxqKsrr+eB6+vSBX74Q7j33uxHkiRJUpuwBKswTJuWvZbb88D1HXccVFRkO0VLkiRJahOWYBWG6dNhiy1g4MDUSdJZf334/vfhzjuz54MlSZIktTpLsNJbtQpmzMiWQpfj88D1nXAC9OjhbLAkSZLURizBSu+xx2DBgvJeCr1ajx5w0knw979n5yZLkiRJalWWYKVXrucDr8kpp0DXrnDuuamTSJIkSSXHEqz0pk/PngUu5+eB66uogAkT4Jpr4LXXUqeRJEmSSoolWGmtWvV/5wPr/5xwAqxcCVdfnTqJJEmSVFIswUpr7lxYuLB8zwdek+22g/32g8sugxhTp5EkSZJKhiVYaa0+H3j//dPmKETHHgtPPQUzZ6ZOIkmSJJWMRktwCKFrCOHBEMKjIYTHQgj/Xz6+VQhhVgjhmRDC5BBCl3x8vfzP8/L3B9b7rjPy8adCCIfWGx+ej80LIUysN97se6jITJ8OW20FW26ZOknh+frXYYMN4PLLUyeRJEmSSkZTZoKXAgfFGHcFhgDDQwh7A+cC58cYBwGLgPH59eOBRTHGbYHz8+sIIewEHAHsDAwHLgwhdAwhdAT+AnwB2Ak4Mr+W5t5DRab++cD6pI02gtGj4frrYfHi1GkkSZKkktBoCY6Z9/M/ds5/InAQcGM+fiVweP77iPzP5O8fHEII+fj1McalMcbngXnAnvnPvBjjczHGZcD1wIj8M829h4pJbS0sWuSmWGszbhy8/z7ceGPj10qSJElqVJOeCc5nbGuAN4G7gWeBt2OMK/JL5gP98t/7AS8D5O+/A/SqP/6xz6xpvFcL7qFi4vnAjdt3Xxg0KNsgS5IkSdI6a1IJjjGujDEOAfqTzdzu2NBl+WtDM7KxFcfXdo+PCCFMCCHMDiHMrqura+AjSmr6dNhmGxgwIHWSwhVCNhtcVQXPPJM6jSRJklT0mrU7dIzxbWA6sDewcQihU/5Wf+DV/Pf5wACA/P0ewML64x/7zJrGF7TgHh/Pe3GMcWiMcWhFRUVz/qpqaytXZs8DOwvcuDFjoEMHuOKK1EkkSZKkoteU3aErQggb5793Az4PPAFMA0bll40Fbs5/vyX/M/n7/4kxxnz8iHxn562AQcCDwEPAoHwn6C5km2fdkn+mufdQsaithbffdlOspth8cxg+PCvBK1emTiNJkiQVtabMBPcFpoUQaskK690xxtuA04FTQwjzyJ7H/Wt+/V+BXvn4qcBEgBjjY8AU4HFgKnBivsx6BXAScBdZuZ6SX0tz76Eisvp8YGeCm2bcOHj1VfjXv1InkSRJkopaKJcJ1KFDh8bZs2enjqHVvvIVePJJePrp1Ek+4RuTZgIw+fh9EiepZ9ky6NcvmzmfMiV1GkmSJKnghBAejjEObey6Zj0TLLWKlSuzjZ6cBW66Ll3g6KPhn/+EBQtSp5EkSZKKliVY7a+mBt55xxLcXOPGwfLlcO21qZNIkiRJRcsSrPbn+cAtM3gwDB0Kf/0rlMljDJIkSVJrswSr/U2bBtttl+16rOYZNy7bWbu6OnUSSZIkqShZgtW+VqyAe+91FriljjgC1lsPLrssdRJJkiSpKFmC1b5qauDddz0fuKV69oSRI+Fvf4MPP0ydRpIkSSo6lmC1r9XnA++/f9ocxWzcOHj7bbj55tRJJEmSpKJjCVb7mj4ddtgB+vZNnaR4HXQQbLGFS6IlSZKkFrAEq/34PHDr6NABjj0W7r4bXnopdRpJkiSpqFiC1X4eeQTee88S3BqOOSY7JunKK1MnkSRJkoqKJVjtx/OBW8/AgXDwwXD55bBqVeo0kiRJUtGwBKv9TJsGO+4Im26aOklpOPZYeP55mDEjdRJJkiSpaFiC1T6WL4f77vNopNY0ciT06OEGWZIkSVIzWILVPh55BN5/36XQralbNzjySPj73+Gdd1KnkSRJkoqCJVjtw/OB28a4cfDBBzB5cuokkiRJUlGwBKt9TJ8OO+8MffqkTlJahg6FXXZxSbQkSZLURJZgtb3VzwO7FLr1hZDNBs+aBY89ljqNJEmSVPAswWp7s2fD4sWW4LZy9NHQqVN2XJIkSZKktbIEq+2tPh/Y54HbRkUFfPnLcNVV2ay7JEmSpDWyBKvtTZuWPbdaUZE6SekaNw7q6uD221MnkSRJkgqaJVhta9kyuP9+zwdua8OHw2abuSRakiRJaoQlWG1r9mxYssTngdtap04wdmw2E/z666nTSJIkSQXLEqy2tfp84M99Lm2OcnDssbByJVx9deokkiRJUsGyBKttTZ8On/oU9O6dOknp23572Hff7MzgGFOnkSRJkgqSJVhtZ/XzwC6Fbj/jxsGTT8IDD6ROIkmSJBUkS7DazoMPwgcfuClWe/r612H99bPZYEmSJEmfYAlW25k+HULweeD2tNFGMHo0XH89LF6cOo0kSZJUcCzBajvTpmXPA2+ySeok5WXcOHj/fbjxxtRJJEmSpIJjCVbbWLoU/vtfl0KnsN9+sO22nhksSZIkNcASrLbx4IPw4YduipVCCNls8IwZMG9e6jSSJElSQbEEq21Mm+bzwCmNGQMdOsAVV6ROIkmSJBUUS7DaxvTpMGQI9OyZOkl56tcPDjkErr4aVq1KnUaSJEkqGJZgtb4PP8yeB3YpdFqVlfDSS3DvvamTSJIkSQXDEqzWN2tWtjGWm2KldfjhsMEG2WywJEmSJMASrLaw+nzgz342dZLytsEG8LWvwQ03wAcfpE4jSZIkFQRLsFrftGmw226w8capk6iyEt59F267LXUSSZIkqSBYgtW6PvwQHnjApdCF4sADYfPNXRItSZIk5SzBal0PPJA9D+ymWIWhY0c46ii4806oq0udRpIkSUrOEqzWNW1adj6tzwMXjspKWLECJk9OnUSSJElKzhKs1jV9Ouy+O/TokTqJVhs8GHbd1SXRkiRJEpZgtaYPPsiWQ7sUuvBUVsKDD8JTT6VOIkmSJCVlCVbrmTULli2zBBeiI4/Mlqlfc03qJJIkSVJSlmC1nocfzl733DNtDn3S5pvDwQdnJXjVqtRpJEmSpGQswWo9NTXQrx9UVKROooZUVsILL8B//5s6iSRJkpSMJVitp7oahgxJnUJr8tWvwvrru0GWJEmSypolWK3jgw/gySdht91SJ9GabLghjBwJU6bAhx+mTiNJkiQlYQlW65g7F1audCa40FVWwttvw+23p04iSZIkJWEJVuuoqclenQkubAcfDH37uiRakiRJZcsSrNZRXQ3du8PAgamTaG06doRvfhPuuAPeeit1GkmSJKndWYLVOmpqsqXQHfyvVMGrrITly2Hy5NRJJEmSpHZnY9G6W7kSHn3UpdDFYtddYfBgl0RLkiSpLFmCte7mzYMlS9wUq5gcfTQ88ED2n50kSZJURizBWnfV1dmrM8HF45vfhBDgmmtSJ5EkSZLalSVY666mBjp3hh13TJ1ETdW/Pxx0UFaCY0ydRpIkSWo3lmCtu+pq2GUX6NIldRI1R2UlPPsszJyZOokkSZLUbhotwSGEASGEaSGEJ0IIj4UQTsnHfx5CeCWEUJP/fLHeZ84IIcwLITwVQji03vjwfGxeCGFivfGtQgizQgjPhBAmhxC65OPr5X+el78/sLF7qJ3FmJVgnwcuPiNHQrdubpAlSZKkstKUmeAVwGkxxh2BvYETQwg75e+dH2Mckv/cAZC/dwSwMzAcuDCE0DGE0BH4C/AFYCfgyHrfc27+XYOARcD4fHw8sCjGuC1wfn7dGu/R4n8KarnXXoO6Op8HLkYbbQRf/Wp2VNLSpanTSJIkSe2i0RIcY3wtxvhI/vt7wBNAv7V8ZARwfYxxaYzxeWAesGf+My/G+FyMcRlwPTAihBCAg4Ab889fCRxe77uuzH+/ETg4v35N91B7q6nJXp0JLk6VlbBoEdxxR+okkiRJUrto1jPB+XLk3YBZ+dBJIYTaEMJlIYSe+Vg/4OV6H5ufj61pvBfwdoxxxcfGP/Jd+fvv5Nev6bvU3lbvDL3rrmlzqGU+/3nYdFOXREuSJKlsNLkEhxA2BP4OfC/G+C5wEbANMAR4Dfj96ksb+HhswXhLvuvjmSeEEGaHEGbX1dU18BGts5oa2GYb6N49dRK1RKdO2XFJt98OCxemTiNJkiS1uSaV4BBCZ7IC/LcY400AMcY3YowrY4yrgEv4v+XI84EB9T7eH3h1LeMLgI1DCJ0+Nv6R78rf7wEsXMt3fUSM8eIY49AY49CKioqm/FXVXNXVPg9c7CorYdkyuOGG1EkkSZKkNteU3aED8FfgiRjjefXG+9a77KvA3Pz3W4Aj8p2dtwIGAQ8CDwGD8p2gu5BtbHVLjDEC04BR+efHAjfX+66x+e+jgP/k16/pHmpP776bHbHj88DFbcgQ2Gknl0RLkiSpLHRq/BL2BSqBOSGEfBckziTb3XkI2TLkF4DjAWKMj4UQpgCPk+0sfWKMcSVACOEk4C6gI3BZjPGx/PtOB64PIZwNVJOVbvLXq0MI88hmgI9o7B5qR48+mr06E1zcQshmg884A557DrbeOnUiSZIkqc00WoJjjPfR8DO4a9xONsb4S+CXDYzf0dDnYozP0cDuzjHGD4GvN+ceakfuDF06jjoKzjwTrrkGfvaz1GkkSZKkNtOs3aGlj6iuhj59oG/fxq9VYRswAA44IFsSHT+xx5wkSZJUMizBarmammwWODS0UEBFp7IS5s2DWbMav1aSJEkqUpZgtcyyZTB3rs8Dl5KvfQ26dnWDLEmSJJU0S7Ba5oknYPlynwcuJd27w+GHw+TJ2b/kkCRJkkqQJVgtU12dvToTXFoqK+Gtt2Dq1NRJJEmSpDZhCVbL1NTA+uvDttumTqLWNGxYttmZS6IlSZJUoizBapnqath1V+jYMXUStaZOneDII+HWW+Htt1OnkSRJklqdJVjNF+P/7Qyt0nP00bB0KdxwQ+okkiRJUquzBKv5nn8e3n3X54FL1R57wA47uCRakiRJJckSrOZbvSmWM8GlKYRsg6x774UXXkidRpIkSWpVlmA1X01N9izwLrukTqK2ctRR2es116TNIUmSJLUyS7Car7o6Wy7brVvqJGorW24J+++fleAYU6eRJEmSWo0lWM1XU+PzwOWgshKeegoeeih1EkmSJKnVWILVPHV18MorPg9cDkaNgq5d3SBLkiRJJcUSrOapqclenQkufT16wIgRcP31sGxZ6jSSJElSq7AEq3ncGbq8VFbCggUwdWrqJJIkSVKrsASreWpqYIstYJNNUidRexg2DPr0cUm0JEmSSoYlWM1TXe0scDnp3BmOPBJuuQUWLUqdRpIkSVpnlmA13eLF2W7BPg9cXsaMyZ4JvuGG1EkkSZKkdWYJVtPNmZOdGetMcHnZbTfYaSe46qrUSSRJkqR1ZglW07kzdHkKIZsNvv9+ePbZ1GkkSZKkdWIJVtNVV0PPntnGWCovRx2VleFrrkmdRJIkSVonlmA1XU1NthQ6hNRJ1N7694eDDsp2iY4xdRpJkiSpxSzBapoVK6C21ueBy1llZbYceubM1EkkSZKkFrMEq2mefho+/NDngcvZyJGw/vqeGSxJkqSiZglW01RXZ6/OBJevjTaCr34VJk+GpUtTp5EkSZJaxBKspqmpgfXWgx12SJ1EKY0ZA4sWwe23p04iSZIktYglWE1TXQ277AKdO6dOopQOPhj69vXMYEmSJBUtS7AaF2M2E+zzwOrYMTsu6Y47YMGC1GkkSZKkZrMEq3Hz58Nbb/k8sDKVlbB8efZssCRJklRkLMFqXE1N9upMsAA+9ansx12iJUmSVIQswWpcdTWEkBUfCbINsmbNgqeeSp1EkiRJahZLsBpXUwODBsGGG6ZOokLxzW9Chw5wzTWpk0iSJEnNYglW46qrfR5YH9W3LxxySLYketWq1GkkSZKkJrMEa+3efhteeMHngfVJlZXw4otw332pk0iSJElNZgnW2q3eFMuZYH3c4YdnS+Q9M1iSJElFxBKstXNnaK3JBhvA174GN9wAH3yQOo0kSZLUJJZgrV11NWy2GWy6aeokKkRjxsC778Itt6ROIkmSJDWJJVhrV1PjLLDW7IADoH9/zwyWJElS0bAEa82WLoXHH/d5YK1Zhw5w9NEwdSq88UbqNJIkSVKjLMFas8cegxUrnAnW2lVWwsqVcP31qZNIkiRJjbIEa82qq7NXZ4K1NjvtBHvs4S7RkiRJKgqWYK1ZTU12BM4226ROokJXWQmPPJKtHpAkSZIKmCVYa1ZdDbvumj33Ka3NkUdCx45ukCVJkqSCZ7tRw1atgkcf9XlgNU2fPjB8OPztb9nzwZIkSVKBsgSrYc8+C++/7/PAaroxY2D+fJg+PXUSSZIkaY0swWpYTU326kywmurLX4bu3V0SLUmSpIJmCVbDqquhUyfYeefUSVQsunWD0aPhxhth8eLUaSRJkqQGWYLVsJqa7Oib9dZLnUTFpLIyK8D//GfqJJIkSVKDLMFqWHW1zwOr+fbbD7bc0jODJUmSVLAswfqk11/PfnweWM3VoUM2G/zvf8Orr6ZOI0mSJH2CJViftHpTLGeC1RKVldkRW9demzqJJEmS9AmWYH2SJVjrYrvtYK+93CVakiRJBckSrE+qroaBA2HjjVMnUbGqrITaWnj00dRJJEmSpI+wBOuTamp8Hljr5hvfgM6dnQ2WJElSwWm0BIcQBoQQpoUQngghPBZCOCUf3ySEcHcI4Zn8tWc+HkIIF4QQ5oUQakMIu9f7rrH59c+EEMbWG98jhDAn/8wFIYTQ0ntoHb3/PjzzjEuhtW5694YvfhH+9jdYsSJ1GkmSJOn/acpM8ArgtBjjjsDewIkhhJ2AicA9McZBwD35nwG+AAzKfyYAF0FWaIGzgL2APYGzVpfa/JoJ9T43PB9v1j3UCh59FGJ0JljrbsyYbJfxe+5JnUSSJEn6fxotwTHG12KMj+S/vwc8AfQDRgBX5pddCRye/z4CuCpmHgA2DiH0BQ4F7o4xLowxLgLuBobn73WPMc6MMUbgqo99V3PuoXXlplhqLYcdBj17wpVXNn6tJEmS1E6a9UxwCGEgsBswC9g0xvgaZEUZ6JNf1g94ud7H5udjaxuf38A4LbiH1lV1NfTqBf37p06iYrfeevDNb8JNN8Fbb6VOI0mSJAHNKMEhhA2BvwPfizG+u7ZLGxiLLRhfa5ymfCaEMCGEMDuEMLuurq6RrxSQzQQPGQKhoX/EUjMdfzwsXepssCRJkgpGk0pwCKEzWQH+W4zxpnz4jdU7SN+zAAAgAElEQVRLkPPXN/Px+cCAeh/vD7zayHj/BsZbco+PiDFeHGMcGmMcWlFR0ZS/anlbvhzmzPF5YLWewYNhn33g4ouzZ80lSZKkxJqyO3QA/go8EWM8r95btwCrd3geC9xcb3xMvoPz3sA7+VLmu4BhIYSe+YZYw4C78vfeCyHsnd9rzMe+qzn30Lp48klYtszngdW6JkyAp56CqqrUSSRJkqQmzQTvC1QCB4UQavKfLwLnAIeEEJ4BDsn/DHAH8BwwD7gEOAEgxrgQ+B/gofznF/kYwHeAS/PPPAvcmY836x5aR9XV2aszwWpNo0dDjx4waVLqJJIkSRKdGrsgxngfDT+DC3BwA9dH4MQ1fNdlwGUNjM8Gdmlg/K3m3kProKYGunaF7bZLnUSlZP31s+OSJk2CBQuyM4QlSZKkRJq1O7RKXHU1fOpT0KnRfzciNc/xx2dL7a+4InUSSZIklTlLsDIxQm0t7Lpr6iQqRTvvDPvu6wZZkiRJSs4SrMyrr8LChZZgtZ3jj4dnnoFp01InkSRJUhmzBCtTW5u9fupTaXOodI0aBT17ZrPBkiRJUiKWYGVWl+DBg9PmUOnq1g3GjoWbboI332z8ekmSJKkNWIKVqa2FAQNg441TJ1EpmzABli93gyxJkiQlYwlWprbWpdBqezvuCJ/9bLYketWq1GkkSZJUhizByo6uefJJS7Dax/HHw7PPwn/+kzqJJEmSypAlWFkBXrHCEqz28bWvQa9eMGlS6iSSJEkqQ5ZguTO02lfXrtkGWf/8J7zxRuo0kiRJKjOWYGUluEsX2G671ElULiZMyFYfXH556iSSJEkqM5ZgZSV4552hU6fUSVQutt8eDjgALrnEDbIkSZLUrizBcmdopTFhAjz3HPz736mTSJIkqYxYgstdXR289hoMHpw6icrNyJHQu7cbZEmSJKldWYLL3Zw52aszwWpv660HxxwDN9+c/YsYSZIkqR1YgsudJVgpTZgAK1fCZZelTiJJkqQyYQkud7W10KcPbLpp6iQqR4MGwUEHuUGWJEmS2o0luNy5KZZSO/54ePFF+Ne/UieRJElSGbAEl7OVK2HuXEuw0jr88Gw1ghtkSZIkqR1YgsvZvHnw4YeWYKXVpQsceyzceiu8+mrqNJIkSSpxluByVlubvXo8klI77rhsZcJf/5o6iSRJkkqcJbic1dZChw6w006pk6jcbbstfP7z2QZZK1emTiNJkqQSZgkuZ3PmwPbbQ9euqZNI2QZZL78MU6emTiJJkqQSZgkuZ+4MrUIyYkR2VNfFF6dOIkmSpBJmCS5X774Lzz9vCVbh6NwZxo2D226D+fNTp5EkSVKJsgSXq7lzs1dLsArJt74FMbpBliRJktqMJbhcrd4Z2hKsQrLVVjBsGFx6KaxYkTqNJEmSSpAluFzV1kKPHjBgQOok0kcdf3y2HPrOO1MnkSRJUgmyBJer2trsfOAQUieRPupLX4LNNoNJk1InkSRJUgmyBJejGLPjkVwKrULUuTOMH5/NBL/0Uuo0kiRJKjGW4HL00kvZ7tCWYBUqN8iSJElSG7EElyM3xVKh23JLGD7cDbIkSZLU6izB5Wh1Cd5ll7Q5pLU5/nh49VW4/fbUSSRJklRCLMHlqLYWtt4aNtoodRJpzQ47DPr1c4MsSZIktSpLcDmqrXUptApfp07ZBllTp8ILL6ROI0mSpBJhCS43H3wATz9tCVZxOO446NAB/vzn1EkkSZJUIizB5eaJJ2DVquyMYKnQDRgAo0fDxRfDO++kTiNJkqQSYAkuN+4MrWJz2mnw3nvZTtGSJEnSOrIEl5vaWujWDbbZJnUSqWn22AMOOAD++EdYvjx1GkmSJBU5S3C5qa3Njkbq2DF1EqnpTjsNXn4ZbrwxdRJJkiQVOUtwOYkRHn3UpdAqPl/8Imy/Pfzud9l/jyVJkqQWsgSXkzfegAULLMEqPh06ZLPBjzwCM2akTiNJkqQiZgkuJ26KpWJWWQkVFfD736dOIkmSpCJmCS4nc+Zkrx6PpGLUtSuceCLcdlt21JckSZLUApbgclJbC5tvDr16pU4itcwJJ2Rl+PzzUyeRJElSkbIEl5PaWpdCq7hVVMDYsXDVVfDmm6nTSJIkqQhZgsvF8uXw+OOWYBW/738fli6Fv/wldRJJkiQVIUtwuXj6aVi2zBKs4rf99vDlL8OFF8IHH6ROI0mSpCJjCS4X7gytUvKDH2THfV11VeokkiRJKjKW4HJRWwudO2ezaFKx++xnYehQOO88WLUqdRpJkiQVEUtwuZgzB3bcEbp0SZ1EWnchwGmnZcv8b7stdRpJkiQVEUtwuait9XxglZZRo2CLLeD3v0+dRJIkSUXEElwOFi2Cl1/2eWCVlk6d4Hvfg6oqeOih1GkkSZJUJCzB5WDOnOzVEqxSM348dO/ubLAkSZKazBJcDtwZWqWqe3eYMAFuvBFefDF1GkmSJBUBS3A5qK2FXr2gb9/USaTW993vZhtl/fGPqZNIkiSpCDRagkMIl4UQ3gwhzK039vMQwishhJr854v13jsjhDAvhPBUCOHQeuPD87F5IYSJ9ca3CiHMCiE8E0KYHELoko+vl/95Xv7+wMbuoTWorc1mgUNInURqfQMGwDe+AZdcAm+/nTqNJEmSClxTZoKvAIY3MH5+jHFI/nMHQAhhJ+AIYOf8MxeGEDqGEDoCfwG+AOwEHJlfC3Bu/l2DgEXA+Hx8PLAoxrgtcH5+3Rrv0by/dhlZtQrmznUptErbaafB++9nRViSJElai0ZLcIyxCljYxO8bAVwfY1waY3wemAfsmf/MizE+F2NcBlwPjAghBOAg4Mb881cCh9f7rivz328EDs6vX9M91JDnn4fFiy3BKm277QYHHZQtiV62LHUaSZIkFbB1eSb4pBBCbb5cumc+1g94ud418/OxNY33At6OMa742PhHvit//538+jV91yeEECaEEGaHEGbX1dW17G9Z7FZviuUZwSp1p50Gr7wCU6akTiJJkqQC1tISfBGwDTAEeA1YfT5JQw+dxhaMt+S7PjkY48UxxqExxqEVFRUNXVL6amuzZ4F33jl1EqltDR8OO+6YHZcUG/x/CZIkSVLLSnCM8Y0Y48oY4yrgEv5vOfJ8YEC9S/sDr65lfAGwcQih08fGP/Jd+fs9yJZlr+m71JDaWhg0CNZfP3USqW116ACnngo1NTBtWuo0kiRJKlAtKsEhhPpn7XwVWL1z9C3AEfnOzlsBg4AHgYeAQflO0F3INra6JcYYgWnAqPzzY4Gb633X2Pz3UcB/8uvXdA81ZPXO0FI5OPpo6NMHfve71EkkSZJUoJpyRNJ1wExg+xDC/BDCeOA3IYQ5IYRa4EDg+wAxxseAKcDjwFTgxHzGeAVwEnAX8AQwJb8W4HTg1BDCPLJnfv+aj/8V6JWPnwpMXNs91vGfQ2lavBiefdYSrPLRtSucdBLceSc8/njqNJIkSSpAIZbJs3NDhw6Ns2fPTh2jfT34IOy1F/zznzBiROo0ReMbk2YCMPn4fRInUYssWABbbAHf/CZcemnqNJIkSWonIYSHY4xDG7tuXXaHVqFbvTO0M8EqJ717wzHHwNVXw+uvp04jSZKkAmMJLmW1tbDhhrDllqmTSO3r+9+H5cvhL39JnUSSJEkFxhJcymprs/OBO/gfs8rMoEHwla/AhRfCkiWp00iSJKmA2I5KVYzuDK3ydtppsHAhXHFF6iSSJEkqIJbgUvXKK7BokSVY5Wu//WDPPeH882GlG8hLkiQpYwkuVW6KpXIXQjYbPG8e3Hpr6jSSJEkqEJbgUjVnTvY6eHDaHFJKI0dmG8P95jfZIwKSJEkqe5bgUlVbm/2P/x49UieR0unUCU4/HWbOhDvuSJ1GkiRJBcASXKrcFEvKHHccbLMNnHEGrFqVOo0kSZISswSXoqVL4cknXQotAXTuDGefnT0icO21qdNIkiQpMUtwKXrySVixwplgabXRo2G33eCnP83+JZEkSZLKliW4FLkztPRRHTrAr38NL7wAF1+cOo0kSZISsgSXotpaWG89GDQodRKpcAwbBgccAP/zP/Dee6nTSJIkKRFLcCmaMwd23jnbGVdSJgQ45xyoq4Pzz0+dRpIkSYlYgkuRO0NLDdtrr+zs4N/+NivDkiRJKjuW4FJTVwevvWYJltbk7LNhyRL41a9SJ5EkSVICluBSM2dO9moJlhq2445w7LFw4YXw4oup00iSJKmdWYJLzeqdoT0jWFqzs87KnhH+2c9SJ5EkSVI7swSXmtpa2HRT6NMndRKpcA0YACefDFdfDXPnpk4jSZKkdmQJLjVuiiU1zcSJsNFGcOaZqZNIkiSpHVmCS8nKlfDYY5ZgqSl69YLTT4dbb4X770+dRpIkSe3EElxK5s2DDz+0BEtNdcopsNlm2axwjKnTSJIkqR1YgkvJ6k2xLMFS02ywQbY51n33wR13pE4jSZKkdmAJLiW1tdCxY3YEjKSmOe442GYbOOOM7JECSZIklTRLcCmprYXtt4f11kudRCoenTvD2WdnZ2xfd13qNJIkSWpjluBS4s7QUsuMHg277QY//SksXZo6jSRJktqQJbhUvPsuvPCCJVhqiQ4d4Ne/zv5vaNKk1GkkSZLUhizBpWLu3OzVEiy1zLBhcOCB2dLo995LnUaSJEltxBJcKqqrs1dLsNQyIcA550BdHZx3Xuo0kiRJaiOW4FJx770wYAD07586iVS89twTRo6E3/0uK8OSJEkqOZbgUhAjzJgBn/tcNpslqeV++UtYsiR7lSRJUsmxBJeCefPg9ddh//1TJ5GK3w47wLHHwkUXZRtlSZIkqaRYgkvBjBnZ6+c+lzaHVCp+/vNsVcVZZ6VOIkmSpFZmCS4FVVWw6aaw3Xapk0iloX9/OPlkuPpqmDMndRpJkiS1IktwKaiq8nlgqbWdcQZ07w4//nHqJJIkSWpFluBi9+KL2Y9LoaXWtckmcPrpcOutcN99qdNIkiSplViCi11VVfZqCZZa33e/C5ttBhMnZruwS5IkqehZgovdjBnQsyfsskvqJFLp2WCDbHOs++/PZoQlSZJU9CzBxa6qCj77Wejgf5RSmxg/HnbaKZsVXrw4dRpJkiStI5tTMXvtNXjmGc8HltpS584waVL27P3Pf546jSRJktaRJbiY+Tyw1D722w++9S04/3yoqUmdRpIkSevAElzMqqpgo41gyJDUSaTSd+650KsXTJgAK1emTiNJkqQWsgQXs6oq2Hdf6NQpdRKp9PXsCX/4Azz0EFx4Yeo0kiRJaiFLcLFasADmznUptNSejjgCDj0UzjwT5s9PnUaSJEktYAkuVvfdl71agqX2E0I2C7xyZbZbtCRJkoqOJbhYzZgBXbvCpz+dOolUXrbeOjs7+B//gJtvTp1GkiRJzWQJLlZVVbDPPtClS+okUvk59VQYPBhOOgneey91GkmSJDWDJbgYvfNOdkyL5wNLaaw+O/iVV+CnP02dRpIkSc1gCS5G998Pq1b5PLCU0j77wLe/DX/6E8yenTqNJEmSmsgSXIyqqrKZqL32Sp1EKm+//jX06ZOdHbxiReo0kiRJagJLcDGaMQP23BPWXz91Eqm89egBF1wA1dXZjLAkSZIKniW42CxenC29dCm0VBhGjYLDDsueDX7ppdRpJEmS1AhLcLF54IFs2aUlWCoMIcCf/wwxwoknZq+SJEkqWJbgYjNjBnToAPvumzqJpNUGDoRf/AJuuw1uuil1GkmSJK2FJbjYVFXB7rvDRhulTiKpvlNOgSFD4OSTs2PMJEmSVJAaLcEhhMtCCG+GEObWG9skhHB3COGZ/LVnPh5CCBeEEOaFEGpDCLvX+8zY/PpnQghj643vEUKYk3/mghBCaOk9St7SpdlyaJdCS4WnUye4+GJ4/XX48Y9Tp5EkSdIaNGUm+Apg+MfGJgL3xBgHAffkfwb4AjAo/5kAXARZoQXOAvYC9gTOWl1q82sm1Pvc8Jbcoyw8+GBWhPffP3USSQ359KfhpJPgwgth1qzUaSRJktSARktwjLEKWPix4RHAlfnvVwKH1xu/KmYeADYOIfQFDgXujjEujDEuAu4GhufvdY8xzowxRuCqj31Xc+5R+qqqsk149tsvdRJJa3L22bD55tnZwcuXp04jSZKkj2npM8GbxhhfA8hf++Tj/YCX6103Px9b2/j8BsZbco/SN2MGDB4Mm2ySOomkNenePdsturYW/vCH1GkkSZL0Ma29MVZoYCy2YLwl9/jkhSFMCCHMDiHMrqura+RrC9zy5fDf//o8sFQMDj8cRoyAs86C559PnUaSJEn1tLQEv7F6CXL++mY+Ph8YUO+6/sCrjYz3b2C8Jff4hBjjxTHGoTHGoRUVFc36Cxac6mpYvNgSLBWLP/0JOnaEE07w7GBJkqQC0tISfAuweofnscDN9cbH5Ds47w28ky9lvgsYFkLomW+INQy4K3/vvRDC3vmu0GM+9l3NuUdpmzEje7UES8VhwIDs+eCpU2HKlNRpJEmSlGvKEUnXATOB7UMI80MI44FzgENCCM8Ah+R/BrgDeA6YB1wCnAAQY1wI/A/wUP7zi3wM4DvApflnngXuzMebdY+SV1UF228Pm26aOomkpjrpJNhjj+wM4UWLUqeRJEkS0KmxC2KMR67hrYMbuDYCJ67hey4DLmtgfDawSwPjbzX3HiVr5Uq4914YPTp1EknN0bFjdnbwpz8NEyfCpEmpE0mSJJW91t4YS21hzhx45x3PB5aK0e67w/e/n5XhW29NnUaSJKnsWYKLQVVV9urzwFJxOvts2G03GDsWXnghdRpJkqSyZgkuBjNmwFZbZRvtSCo+XbvCDTdkjzaMHg3LlqVOJEmSVLYswYUuxmwm2Flgqbhtsw1cfjk89BD88Iep00iSJJUtS3Che+IJWLDAEiyVgpEjs52iL7gAbrwxdRpJkqSyZAkudKufB3ZTLKk0/OY3sOeeMH48zJuXOo0kSVLZsQQXuqoq2Hxz2Hrr1EkktYYuXWDKlOz4pK9/HT78MHUiSZKksmIJLmQxZptife5zEELqNJJay5ZbwlVXQU0NfO97qdNIkiSVFUtwIXvuOXj1VZdCS6XoS1+CH/0IJk2Ca69NnUaSJKlsWIILmecDS6Xt7LNhv/1gwgR48snUaSRJksqCJbiQzZgBvXvDjjumTiKpLXTuDNdfD926wahRsGRJ6kSSJEklzxJcyFafD+zzwFLp6tcP/vY3ePxxOPHE1GkkSZJKniW4UL38Mjz/vEuhpXIwbBj85CdwxRVw+eWp00iSJJU0S3Ch8nxgqbycdRYceGA2GzxnTuo0kiRJJcsSXKiqqqBHDxg8OHUSSe2hY8dsl+gePbLzg997L3UiSZKkkmQJLlQzZmS7xnbsmDqJpPay2WZw3XXwzDPw7W9nZ4VLkiSpVVmCC9Ebb8BTT7kUWipHBxwAv/hFNit88cWp00iSJJUcS3Ah8nxgqbydcQYceih897vwyCOp00iSJJUUS3AhqqqCDTaA3XdPnURSCh06wDXXQEVF9nzwO++kTiRJklQyLMGFqKoKPvMZ6Nw5dRJJqfTuDZMnw4svwvjxPh8sSZLUSizBhWbhwux4FJdCS9p3XzjnHPj73+FPf0qdRpIkqSRYggvNffdlMz5uiiUJ4LTT4Mtfhh/8AB54IHUaSZKkomcJLjQzZsB668GnP506iaRCEAJceSX07w9f+Qo8/XTqRJIkSUXNElxoqqpgr72ga9fUSSQVip494a67st8POQReeSVtHkmSpCJmCS4k772XHYfiUmhJHzdoEEydCosWwbBh2f4BkiRJajZLcCG5/35YtcpNsSQ1bPfd4eabYd48OOwwWLw4dSJJkqSiYwkuJFVV0KkT7LNP6iSSCtWBB8L118ODD8KoUbBsWepEkiRJRcUSXEiqqmDoUNhgg9RJJBWyr34VJk3Klkcfc0y2gkSSJElNYgkuFEuWZDM7LoWW1BTHHQe//jVcdx1873vZ0WqSJElqVKfUAZSbNQuWL3dTLElNd/rpUFcH550HFRXw05+mTiRJklTwLMGFYsaM7DzQffdNnURSsQgBfvtbWLAAfvYz6N0bvvOd1KkkSZIKmiW4UFRVwZAh0KNH6iSSikmHDnDppdmRSSeeCL16wejRqVNJkiQVLJ8JLgTLlsHMmS6FltQynTvDlCnZSpKjj4a7706dSJIkqWBZggtBx47Z/2g9/vjUSSQVq27d4NZbYccds92jH3wwdSJJkqSCZAkuBB07wn77wQ47pE4iqZhtvHF2bFKfPvDFL8ITT6ROJEmSVHAswZJUSvr2hX/9Czp1gmHD4OWXUyeSJEkqKJZgSSo1226bzQi/+25WhBcsSJ1IkiSpYFiCJakUDRkCt9wCzz8Phx0G77+fOpEkSVJBsARLUqnaf3+YPBlmz4aRI2Hp0tSJJEmSkrMES1IpGzEiO0f47rthzBhYsSJ1IkmSpKQ6pQ4gSWpjxx4Lb70FP/whLFkC118PG2yQOpUkSVISzgRLUjn4wQ/goovgjjvgoIOgri51IkmSpCQswZJULr79bbjpJqithX33heeeS51IkiSp3VmCJamcjBgB99yTLY/eZx94+OHUiSRJktqVJViSys1nPgP33w/dumU7SE+dmjqRJElSu7EES1I52mEHmDkTBg2CL38ZrrwydSJJkqR2YQmWpHLVty/MmAEHHADHHAO/+hXEmDqVJElSm7IES1I5694dbr8djjoKfvxjOPFEWLkydSpJkqQ24znBklTuunSBq66C/v3h3HPhtdfg2muzZ4YlSZJKjDPB+v/bu/Mgu8oyj+Pfp7uTkJCFBLKRlSWsYTUsBQwEZAkIgiLKOAhaIuAwiKUUi6WyI6KoY4EYNkUqAxOV1UECQca4AoEECJsJECAEEkzYAgGSzjt/vPdOL+lOOkl3n9t9v5+qt8457z333ifw6uWXc877ShLU1MDll8NPfwp33gkHH5xnkJYkSepmDMGSpAZnnAFTp+alk/bdF+bPL7oiSZKkdmUIliQ19ZnPwH33waJFeS3h2bOLrkiSJKndGIIlSavbf3/485+hri7vT59edEWSJEntwhAsSWrZjjvmtYTHjIHDD4cpU4quSJIkaYMZgiVJrRs5Ev70J9hvPzjhBLj4YpdQkiRJXdoGheCImB8RT0bE7IiYWeobFBH3R8Tc0nZgqT8i4qcRMS8inoiI3Rt9zkml8+dGxEmN+j9W+vx5pffGmr5DktQBNtkE7r03ryX83e/CYYfB668XXZUkSdJ6aY8rwQemlHZNKU0oHZ8LPJBSGgc8UDoGOBwYV2qnANdADrTA+cBewJ7A+Y1C7TWlc8vvm7SW75AkdYReveDmm+H66+Gvf4VddoFp04quSpIkaZ11xO3QRwM3lfZvAo5p1P+rlP0d2CQihgOHAfenlJamlN4E7gcmlV7rn1L6W0opAb9q9lktfYckqaNEwJe/DDNnwpAhMGkSnHMOrFhRdGWSJElttqEhOAH3RcSjEXFKqW9oSuk1gNJ2SKl/BPBKo/cuKPWtqX9BC/1r+g5JUkfbYQd4+GE49VS44gr4l3+BF18suipJkqQ22dAQvG9KaXfyrc6nR8T+azg3WuhL69HfZhFxSkTMjIiZb7zxxrq8VZK0Jr17w89/DlOnwrPPwm67wa9/XXRVkiRJa7VBITiltLC0XQzcTn6md1HpVmZK28Wl0xcAoxq9fSSwcC39I1voZw3f0by+a1NKE1JKEwYPHry+f0xJUmuOOw5mzYLttoPPfhZOOw2WLy+6KkmSpFatdwiOiI0jol95HzgUmAPcBZRneD4JuLO0fxdwYmmW6L2Bt0u3Mk8DDo2IgaUJsQ4FppVeezci9i7NCn1is89q6TskSZ1tiy3yMkrnnAOTJ8Oee8LTTxddlSRJUos25ErwUODPEfE48DDwPymle4HLgUMiYi5wSOkY4B7gBWAecB3w7wAppaXAxcAjpXZRqQ/gq8D1pfc8D/y+1N/ad0iSitCjB1x+eV5KafFimDABbrgB0jo9xSJJktTh6tb3jSmlF4BdWuhfAny8hf4EnN7KZ90I3NhC/0xgfFu/Q5JUsMMOg8cfhy98AU4+GaZPz88ODxhQdGWSJElAxyyRJEmqZsOG5TWEL7ssT5a1++7wyCNFVyVJkgQYgiVJHaGmBs47D2bMgJUrYZ994MorYdWqoiuTJElVzhAsSeo4++wDs2fDUUfBWWfBkUfCK6+s/X2SJEkdxBAsSepYAwfCb38LP/sZPPhgXk7pkktcSkmSJBXCECxJ6ngR8NWvwrPPwhFHwHe+AzvsALfd5gzSkiSpUxmCJUmdZ8yYPFnWH/4A/frBscfCwQfDnDlFVyZJkqqEIViS1PkOPBAeewyuugpmzYJdd4UzzoClS9f+XkmSpA1gCJYkFaOuDk4/HebOhVNPzc8Mb7NNXle4vr7o6iRJUjdlCJYkFWvTTeHqq/MV4fHj87PDH/tYXl5JkiSpnRmCJUmVYeed8+zRU6fCm2/CAQfA5z4HL79cdGWSJKkbMQRLkipHBBx3HDzzDJx/Ptx1V15S6aKLXFJJkiS1C0OwJKny9OkDF1yQl1Q68sgciLffHn7zG5dUkiRJG8QQLEmqXGPG5NujH3wQBgzIV4n33x+mTTMMS5Kk9WIIliRVvokT4dFH8wzSL74IkybBHnvAbbfBqlVFVydJkroQQ7AkqWuoq8szRz//PFx3Hbz1Fhx7bJ5R+uabYcWKoiuUJEldgCFYktS19OoFJ5+cnxe+5ZYcjk88sWGN4Q8+KLpCSZJUwQzBkqSuqa4Ojj8eZs/Os0gPHZqvFG+5JVx5JSxbVnSFkiSpAhmCJUldW00NHHUU/O1v8MADsMMOcNZZeVKtCy+EpUuLrlCSJFUQQ7AkqXuIgIMOgunTcyDeb7+8zM4fyHQAAA3eSURBVNKYMXD22fD660VXKEmSKoAhWJLU/ey9N9x5JzzxBHzyk/n26LFj4fTTYf78oquTJEkFMgRLkrqvnXaCKVPguefy5FnXXQdbb51nlb73XqivL7pCSZLUyQzBkqTub+ut4dpr4YUX4BvfgBkz4PDD8yRaF14IL79cdIWSJKmTGIIlSdVj5Ei44gp49VWYOhW22y6H4LFj4Ygj4PbbXW9YkqRuzhAsSao+PXvCccfBtGn56vC3v52fH/70p2HUKDj3XJg7t+gqJUlSBzAES5Kq29ixcNFFecKsu+/Ok2r98IewzTZw4IH5meIPPii6SkmS1E4MwZIkAdTVwZFHwh135GeEL7ssb084ATbfHL72NXjyyaKrlCRJG8gQLElSc5tvDuedl2+Jnj4dDjsMJk+GnXfOV4onT4Y33ii6SkmStB4MwZIktaamBj7+cbjlljyZ1o9/DO++C6edBsOG5dulr7oKFi4sulJJktRGhmBJktpis83g61+HOXNg1iz41rdg0SI44wwYMQL23Rd+9KP8bLEkSapYhmBJktZFBOy6K1x8MTz9dG4XXwzvvw/f/CZssQVMmADf+x784x9FVytJkpoxBEuStCG23z4vsTRrFsybB9//PtTW5ivF224LO+0EF1yQJ9VKqehqJUmqeoZgSZLay1Zbwdlnw0MP5Zmlf/ITGDgwL8G08845FJ93HsycaSCWJKkghmBJkjrCqFFw5pkwY0aeOOuaa2DMGPjBD2CPPWD48Lz80i9/CQsWFF2tJElVo67oAiRJ6vaGDcszSp92GixZAnffDffdl9uUKfmc7baDgw+GQw6BAw6AAQOKrVmSpG7KECxJUmfadFP44hdzW7UqPys8fXpuN9yQl1yqrYU992wIxXvtBT17Fl25JEndgrdDS5JUlJoa2GWXPKv0738Pb74JDz4I556bA/Kll8L++8OgQfCJT+R1iufM8XliSZI2gFeCJUmqFL16wcSJuV1yCbz1Vg7F5SvF99yTzxs6FA46KF8h3mMP2G036N27yMolSeoyDMGSJFWqTTaBT30qN8gzTj/wANx/P/zxj3DLLbm/rg7Gj8+3UO+xR2477pj7JUlSE/46SpLUVYweDV/6Um4Ar74KjzzS0KZOhWuvza/17g27754DcTkcb7UVRBRXvyRJFcAQLElSVzViRG7HHJOPU4J58xpC8cMPw+TJeb1iyGsWl68Ul2+jHjXKYCxJqiqGYEmSuosIGDcut89/PvetXAlPPZUDcTkcX3451Nfn1/v1y7dOjx+ft+X9YcMMx5KkbskQLElSd1ZXl2eg3mUX+MpXct/778Ps2fDEEzkgz5kDd9wB11/f8L6BAxuCceOAPHhwMX8OSZLaiSFYkqRq06cP7LNPbo0tXtwQisvbW2/Ns1SXDRnSNBhvu21uXjmWJHURhmBJkpQNGZLbgQc29KUEr73WNBg/9RT84hewbFnDef37wzbbNITichs3LoduSZIqhCFYkiS1LgI23zy3Qw9t6F+1Cl55BZ57rmmbMQOmTGn6GaNHrx6Ot90WRo6EmprO/fNIkqqeIViSJK27mhoYMya3xuEY4L338izVzQPyTTfBu+82nNenD2y5ZcPnjB3bsD9mDAwd6i3WkqR2ZwiWJEnta+ONGybjaiwleP31psH4hRdg/nz4y1+aPnsM0KtXvorcWkgeMSJP/CVJ0jrwl0OSJHWOCBg+PLeJE1d//Z134KWXGtr8+Q37v/sdLFrU9Pza2hyER45suGW73EaMaNjv188rypKk/2cIliRJlaF/f9hpp9xasnw5vPxy06D80kuwcGGesOu++3KQbm7jjZuG4uZBefjwPCFY376GZUmqAoZgSZLUNfTu3TCpVmuWLcuhuHF79dWG/b//PR9/+OHq791oo4YZsstt8ODV+8r9vXp13J9VktRhDMGSJKn76Ns3L9W0zTatn5NSfv64HI5few3eeCOvk1xuixbBk0/m/ZYCM8CAAQ2heNNNcxs0qOm2eZ/LRUlS4QzBkiSpukTAwIG5jR+/5nNTyjNaNw7I5VYOzosW5eeXH3sMlizJt223ZqONWg7LgwblUF1u/fs3PS731da26z8KSapGhmBJkqTWROTw2b8/bL11296zfDksXZoD8ZIlLe+Xt888k/eXLoUVK9b+2X37rh6Oy61fv4bWt++aj/v0cY1mSVXLECxJktSeevfOk26NGNH296QEH3wAb7/d0N55p+lxS/1LluRlpt5+O1+xfv/9tn1fRJ4wrHlI3njjDWt9+kCPHk4wJqmiGYIlSZKKFpHDc+/eMGzY+n9OfT28914OxOW2bFnrx81fW7Qov79xa+2Z6Lb8WTbaqOm2LfuNt+uy37On4VtSm3TpEBwRk4D/BGqB61NKlxdckiRJUnFqaxtu324vK1fmK8zNw3FL7f338xXt5ctza2n/7bfh9ddXf335cli1asNqLYfinj3z7N3l1vi4ra+V93v2bLrflm2PHrnV1TVs6+oM6VKF6LIhOCJqgauBQ4AFwCMRcVdK6eliK5MkSepG6uraP1i3ZsWKHIzL4Xh99z/6KF/BLrfGx+++C//85+r9jY83NIy3prZ29XDc2rZ5q61d97617bf19dra1vfX9vraWk2NfzmgTtdlQzCwJzAvpfQCQETcChwNGIIlSZK6ovIV1H79iq2jvr5pKG5pu6bXPvwwX0EvtxUr1n1bX9/0Mz76aPW+lStb72v8GfX1HRfs20NNzdqD8oa08mc0/rzW+tZ0bvP9NfU1749ofbum1xqf0/zctrbm39O8prYcl/e32y4/ftDFdeUQPAJ4pdHxAmCvgmqRJElSd1Fbmyf56k7rOqfUNDC3tN88UNfXN32t+XltfX192qpVTffXtzX/rPJfJjTub/5d61JDa691V3PmwI47Fl3FBuvKIbil+yZSkxMiTgFOARg9enRn1CRJkiRVnoiG253V8ZoH5ZRyW7Wq5e2aXmt8TvNz29qaf0+5rctxSjBqVNH/ZNtFV/5fwQKg8b+FkcDCxieklK4FrgWYMGFCk4AsSZIkSR2ifBuxKlJX/jfzCDAuIraIiJ7A8cBdBdckSZIkSapgXfZKcEppZUT8BzCNvETSjSmlpwouS5IkSZJUwbpsCAZIKd0D3FN0HZIkSZKkrqEr3w4tSZIkSdI6MQRLkiRJkqqGIViSJEmSVDUMwZIkSZKkqmEIliRJkiRVDUOwJEmSJKlqGIIlSZIkSVXDECxJkiRJqhqGYEmSJElS1TAES5IkSZKqhiFYkiRJklQ1DMGSJEmSpKphCJYkSZIkVQ1DsCRJkiSpahiCJUmSJElVwxAsSZIkSaoahmBJkiRJUtWIlFLRNXSKiHgDeKnoOtZiM+CfRRchNeKYVCVyXKoSOS5VaRyTqkQdPS7HpJQGr+2kqgnBXUFEzEwpTSi6DqnMMalK5LhUJXJcqtI4JlWJKmVceju0JEmSJKlqGIIlSZIkSVXDEFxZri26AKkZx6QqkeNSlchxqUrjmFQlqohx6TPBkiRJkqSq4ZVgSZIkSVLVMARXgIiYFBHPRcS8iDi36HpUnSLixohYHBFzGvUNioj7I2JuaTuwyBpVXSJiVEQ8GBHPRMRTEXFmqd9xqcJExEYR8XBEPF4alxeW+reIiIdK4/K/I6Jn0bWq+kREbUTMiojflY4dlypMRMyPiCcjYnZEzCz1VcRvuCG4YBFRC1wNHA7sAPxrROxQbFWqUr8EJjXrOxd4IKU0DnigdCx1lpXAN1NK2wN7A6eX/v/RcakifQgclFLaBdgVmBQRewPfB35cGpdvAl8usEZVrzOBZxodOy5VtANTSrs2WhapIn7DDcHF2xOYl1J6IaX0EXArcHTBNakKpZRmAEubdR8N3FTavwk4plOLUlVLKb2WUnqstP8u+T/sRuC4VIFStqx02KPUEnAQ8JtSv+NSnS4iRgKfAK4vHQeOS1WeivgNNwQXbwTwSqPjBaU+qRIMTSm9BjmQAEMKrkdVKiLGArsBD+G4VMFKt5zOBhYD9wPPA2+llFaWTvG3XEX4CXA2sKp0vCmOSxUrAfdFxKMRcUqpryJ+w+uK+FI1ES30OWW3JJVERF/gt8DXU0rv5IsbUnFSSvXArhGxCXA7sH1Lp3VuVapmEXEksDil9GhETCx3t3Cq41Kdad+U0sKIGALcHxHPFl1QmVeCi7cAGNXoeCSwsKBapOYWRcRwgNJ2ccH1qMpERA9yAJ6SUrqt1O24VEVIKb0F/C/5mfVNIqJ8ccHfcnW2fYFPRsR88qN1B5GvDDsuVZiU0sLSdjH5Lwz3pEJ+ww3BxXsEGFeava8ncDxwV8E1SWV3ASeV9k8C7iywFlWZ0vNsNwDPpJR+1Oglx6UKExGDS1eAiYjewMHk59UfBD5TOs1xqU6VUjovpTQypTSW/N+Sf0gp/RuOSxUkIjaOiH7lfeBQYA4V8hseKXlXRNEi4gjy39bVAjemlC4tuCRVoYi4BZgIbAYsAs4H7gCmAqOBl4HjUkrNJ8+SOkRE7Af8CXiShmfcvkV+LthxqUJExM7kyVxqyRcTpqaULoqILclX4AYBs4ATUkofFlepqlXpduizUkpHOi5VlNLYu710WAf8V0rp0ojYlAr4DTcES5IkSZKqhrdDS5IkSZKqhiFYkiRJklQ1DMGSJEmSpKphCJYkSZIkVQ1DsCRJkiSpahiCJUmSJElVwxAsSZIkSaoahmBJkiRJUtX4PynJ6XrIyxMWAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1152x648 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "dps = list(sample_dp.keys())\n",
    "dps.sort()\n",
    "dp_dist = [sample_dp[x] for x in dps]\n",
    "fig, ax = plt.subplots(figsize=(16, 9))\n",
    "ax.plot(dp_dist[:50], 'r')\n",
    "ax.axvline(dp_dist.index(max(dp_dist)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
